/*
 * GeoDatabase.cpp
 *
 *  Created on: 13 Jul 2011
 *      Author: allan
 */

#include "GeoDatabase.h"

// C++ includes
#include <algorithm>
#include <cstdlib>
#include <fstream>
using namespace std;

// Boost includes
#include <boost/algorithm/string.hpp>
#include <boost/foreach.hpp>
using namespace boost;

// GeoReact includes
#include <Utils/Algorithms.h>

string 
GetLine(ifstream& file)
{
	string line; 
	
	std::getline(file, line); 
	
	boost::trim(line);
	
	return line;
}

double
StringToFloat(string num)
{
	return atof(num.c_str());
}

vector<string> 
SplitString(string word, string tokens)
{
	vector<string> words;
	
	boost::split(words, word, boost::is_any_of(tokens));
	
	boost::erase_all(words, vector<string>(1, ""));
	
	BOOST_FOREACH(string& w, words)
		boost::trim(w);
	
	return words;
}

vector<pair<double,string>> 
ParseElementalComposition(string value)
{
	vector<pair<double,string>> elemental_composition;
	
	vector<string> words = SplitString(value, ",");
	
	BOOST_FOREACH(string word, words)
	{
		vector<string> pair = SplitString(word, "()");
		
		double num_atoms = StringToFloat(pair[1]);
		string element   = pair[0];
		
		elemental_composition.push_back({num_atoms, element});
	}
	
	return elemental_composition;
}

double
ParseQuantityWithUnit(string value)
{
	vector<string> words = SplitString(value, " ");
	
	double quantity = StringToFloat(words[0]);
	//string unit     = words[1];
	
	return quantity;
}

vector<pair<double,double>> 
ParseEquilibriumConstants(string value)
{
	vector<pair<double,double>> logK;
	
	vector<string> words = SplitString(value, ",");
	
	BOOST_FOREACH(string word, words)
	{
		vector<string> pair = SplitString(word, "()");
		
		double eq_constant = StringToFloat(pair[0]);
		double temperature = StringToFloat(pair[1]);
		
		logK.push_back({eq_constant, temperature});
	}
	
	return logK;
}

vector<pair<double,string>> 
ParseReactionProducts(string value)
{
	vector<pair<double,string>> products;
	
	vector<string> words = SplitString(value, ",");
	
	BOOST_FOREACH(string word, words)
	{
		vector<string> pair = SplitString(word, " ");
		
		double stoichiometry = StringToFloat(pair[0]);
		string product       = pair[1];
		
		products.push_back({stoichiometry, product});
	}
	
	return products;
}

void ReadBasisSpecies(string basis_filename, set<string>& basis)
{
	// Open the file with the basis species 
	ifstream file(basis_filename);
	
	assert(file.is_open());
	
	while(!file.eof())
	{
		string line = GetLine(file);
		
		if(line == "Basis:" || line == "" || line[0] == '#')
			continue;
		
		basis.insert(line);		
	}
}

void ReadSpecies(string species_filename, map<string, SpeciesData>& species)
{
	// Open the file with the species 
	ifstream file(species_filename);
	
	assert(file.is_open());
	
	while(!file.eof())
	{
		string line = GetLine(file);
		
		if(line == "" || line[0] == '#')
			continue;
		
		if(line == "Species:")
		{
			SpeciesData species_data;
			
			while(!file.eof())
			{
				string line = GetLine(file);
				
				if(line == "")
					break;
				
				if(line[0] == '#')
					continue;
				
				vector<string> words = SplitString(line, "=");
				
				string attribute = words[0];
				string value     = words[1];
				
				if(attribute == "name")
					species_data.name = value;
				
				if(attribute == "phase")
					species_data.phase = value;
				
				if(attribute == "elemental_composition")
					species_data.elemental_composition = ParseElementalComposition(value);
				
				if(attribute == "molecular_weight")
					species_data.molecular_weight = ParseQuantityWithUnit(value);
				
				if(attribute == "volumetric_weigth")
					species_data.volumetric_weigth = ParseQuantityWithUnit(value);
				
				if(attribute == "radius")
					species_data.radius = ParseQuantityWithUnit(value);
			}
			
			species[species_data.name] = species_data;
		}
	}
}

void ReadReactions(string reactions_filename, map<string, ReactionData>& reactions)
{
	// Open the file with the reactions
	ifstream file(reactions_filename);
	
	assert(file.is_open());
	
	while(!file.eof())
	{
		string line = GetLine(file);
		
		if(line == "" || line[0] == '#')
			continue;
		
		if(line == "Reaction:")
		{
			ReactionData reaction_data;
			
			while(!file.eof())
			{
				string line = GetLine(file);
				
				if(line == "")
					break;
				
				if(line[0] == '#')
					continue;
				
				vector<string> words = SplitString(line, "=");
				
				string attribute = words[0];
				string value     = words[1];
				
				if(attribute == "reactant")
					reaction_data.reactant = value;
				
				if(attribute == "products")
					reaction_data.products = ParseReactionProducts(value);
				
				if(attribute == "log(K)")
					reaction_data.logK = ParseEquilibriumConstants(value);
			}
			
			reactions[reaction_data.reactant] = reaction_data;
		}
	}
}

void ReadComplexes(string complexes_filename, map<string, ComplexData>& complexes)
{
	// Open the file with the reactions
	ifstream file(complexes_filename);
	
	assert(file.is_open());
	
	while(!file.eof())
	{
		string line = GetLine(file);
		
		if(line == "" || line[0] == '#')
			continue;
		
		if(line == "Complex:")
		{
			ComplexData complex_data;
			
			while(!file.eof())
			{
				string line = GetLine(file);
				
				if(line == "")
					break;
				
				if(line[0] == '#')
					continue;
				
				vector<string> words = SplitString(line, "=");
				
				string attribute = words[0];
				string value     = words[1];
				
				if(attribute == "name")
					complex_data.name = value;
				
				if(attribute == "composition")
					complex_data.composition = ParseReactionProducts(value);
			}
			
			complexes[complex_data.name] = complex_data;
		}
	}
}

// Initialize the static data members of class GeoDatabase
std::shared_ptr<GeoDatabase> 
GeoDatabase::instance = std::shared_ptr<GeoDatabase>();
//GeoDatabase::instance = 0;

string 
GeoDatabase::path = "Databases/";

GeoDatabase::GeoDatabase()
{}

GeoDatabase::GeoDatabase(const string& path)
{
	string basis_filename = path;
	       basis_filename.append("geo-basis.dtb");
	
	string species_filename = path;
	       species_filename.append("geo-species.dtb");
	
	string reactions_filename = path;
	       reactions_filename.append("geo-reactions.dtb");
	
	string complexes_filename = path;
	       complexes_filename.append("geo-complexes.dtb");
	
	ReadBasisSpecies(basis_filename, basis);
	
	ReadSpecies(species_filename, species);
	
	ReadReactions(reactions_filename, reactions);
	
	ReadComplexes(complexes_filename, complexes);
}

void
GeoDatabase::SetDatabasePath(const string& path_)
{
	path = path_;
	
	instance = std::shared_ptr<GeoDatabase>(new GeoDatabase(path_));
}

const string
GeoDatabase::GetDatabasePath()
{
	return path;
}

const GeoDatabase& 
GeoDatabase::GetInstance()
{
	if(!instance)
		SetDatabasePath(path);
	
	return *instance;
}

const GeoDatabase::ReactionSet&
GeoDatabase::GetReactionSet() const
{
	return reactions;
}

const GeoDatabase::SpeciesSet&
GeoDatabase::GetSpeciesSet() const
{
	return species;
}

const GeoDatabase::ComplexSet&
GeoDatabase::GetComplexSet() const
{
	return complexes;
}

const GeoDatabase::BasisSet&
GeoDatabase::GetBasisSet() const
{
	return basis;
}

const boost::optional<ReactionData>
GeoDatabase::FindReaction(const string& reactant) const
{
	auto iter = reactions.find(reactant);
	
	if(iter != reactions.end())
		return boost::optional<ReactionData>(iter->second);
	else
		return boost::optional<ReactionData>();
}

const boost::optional<SpeciesData>
GeoDatabase::FindSpecies(const string& species_name) const
{
	auto iter = species.find(species_name);
	
	if(iter != species.end())
		return boost::optional<SpeciesData>(iter->second);
	else
		return boost::optional<SpeciesData>();
}

const boost::optional<ComplexData>
GeoDatabase::FindComplex(const string& complex) const
{
	auto iter = complexes.find(complex);
	
	if(iter != complexes.end())
		return boost::optional<ComplexData>(iter->second);
	else
		return boost::optional<ComplexData>();
}

const bool
GeoDatabase::IsBasis(const string& species)
{
	const BasisSet& basis = GetInstance().basis;
	
	return basis.find(species) != basis.end();
}

const bool
GeoDatabase::IsSecondary(const string& species)
{
	return !IsBasis(species);
}

const bool
GeoDatabase::IsComplex(const string& species)
{
	const ComplexSet& complexes = GetInstance().complexes;
	
	return complexes.find(species) != complexes.end();
}

const bool
GeoDatabase::IsIonic(const string& species)
{
	return species.find("+]") != species.npos || 
	       species.find("-]") != species.npos;
}

const bool
IsPhaseTypeHelper(const string& species, const string& phase_type)
{
	auto opt_species = GeoDatabase::GetInstance().FindSpecies(species);
	
	if(opt_species)
		return opt_species->phase == phase_type;
	else
		return false;
}

const bool
GeoDatabase::IsAqueous(const string& species)
{
	return IsPhaseTypeHelper(species, "aqueous");
}

const bool
GeoDatabase::IsGaseous(const string& species)
{
	return IsPhaseTypeHelper(species, "gaseous");
}

const bool
GeoDatabase::IsMineral(const string& species)
{
	return IsPhaseTypeHelper(species, "mineral");
}

const double
GeoDatabase::ElectricalCharge(const string& species)
{
	const size_t i1 = species.find("[");
	const size_t i2 = species.find("]");
	
	if(i1 == string::npos || i2 == string::npos)
		return 0;
	
	string numeral = (i2-i1 > 2) ? species.substr(i1+1, i2-i1-1) : "1";
	string sign = species.substr(i2-1, 1);
	
	string charge = sign.append(numeral);
	
	return atof(charge.c_str());
}

const double
GeoDatabase::ElementAtoms(const string& species, const string& element)
{
	auto opt_species = GeoDatabase::GetInstance().FindSpecies(species);
	
	if(!opt_species)
		return 0;
	
	BOOST_FOREACH(auto pair, opt_species->elemental_composition)
		if(pair.second == element)
			return pair.first;
	
	return 0;
}

const double
GeoDatabase::MolecularWeight(const string& species)
{
	auto opt_species = GeoDatabase::GetInstance().FindSpecies(species);
	
	return (opt_species) ? opt_species->molecular_weight : 0.0; 
}

const double
GeoDatabase::VolumetricWeight(const string& species)
{
	auto opt_species = GeoDatabase::GetInstance().FindSpecies(species);
	
	return (opt_species) ? opt_species->volumetric_weigth: 0.0; 
}

const vector<string>
SpeciateSpecies(const vector<string>& components)
{
	const vector<string> basis_species = SpeciateBasisSpecies(components);
	
	const vector<string> secondary_species = SpeciateSecondarySpecies(components);
	
	return Union(basis_species, secondary_species);
}

const vector<string>
SpeciateAqueousSpecies(const vector<string>& components)
{
	return Filter(SpeciateSpecies(components), GeoDatabase::IsAqueous);
}

const vector<string>
SpeciateGaseousSpecies(const vector<string>& components)
{
	return Filter(SpeciateSpecies(components), GeoDatabase::IsGaseous);
}

const vector<string>
SpeciateMineralSpecies(const vector<string>& components)
{
	return Filter(SpeciateSpecies(components), GeoDatabase::IsMineral);
}

const vector<string>
SpeciateBasisSpecies(const vector<string>& components)
{
	set<string> basis_species;
	
	// Make sure the hydrogen ion H[+] is present in the set of primary species in case H2O is listed in (components)
	if(IsContained("H2O", components))
		basis_species.insert("H[+]");
	
	BOOST_FOREACH(const string& component, components)
	{
		if(GeoDatabase::IsBasis(component))
			basis_species.insert(component);
		else
		{
			auto opt_reaction_data = GeoDatabase::GetInstance().FindReaction(component); assert(opt_reaction_data);
			
			// Insert the products of the reaction in the set of basis species 
			BOOST_FOREACH(auto pair, opt_reaction_data->products)
				if(GeoDatabase::IsBasis(pair.second))
					basis_species.insert(pair.second);
		}
	}
	
	return vector<string>(basis_species.begin(), basis_species.end());
}

const vector<string>
SpeciateSecondarySpecies(const vector<string>& components)
{
	const vector<string> basis_species = SpeciateBasisSpecies(components);
	
	vector<string> secondary_species;
	
	BOOST_FOREACH(auto reaction, GeoDatabase::GetInstance().GetReactionSet())
	{
		const string& reactant = reaction.first;
		
		const ReactionData& reaction_data = reaction.second;
		
		bool all_products_are_basis = true;
		
		BOOST_FOREACH(auto product, reaction_data.products)
			if(!IsContained(product.second, basis_species))
				all_products_are_basis = false;
		
		if(all_products_are_basis)
			secondary_species.push_back(reactant);
	}
	
	return secondary_species;
}
